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2 . ABSTRACT 

■ We use semi-analytic modeling on top of the Millennium simulation to study the joint forma- 
ts \ tion of galaxies and their embedded supermassive black holes. Our goal is to test scenarios 

in which black hole accretion and quasar activity are triggered by galaxy mergers, and to 
constrain different models for the lightcurves associated with individual quasar events. In the 
present work we focus on studying the spatial distribution of simulated quasars. At all lumi- 
nosities, we find that the simulated quasar two-point correlation function is fit well by a single 
' power-law in the range 0.5 < r < 20/i^'Mpc, but its normalization is a strong function of 

, redshift. When we select only quasars with luminosities within the range typically accessi- 

■ ble by today's quasar surveys, their clustering strength depends only weakly on luminosity, 

in agreement with observations. This holds independently of the assumed lightcurve model, 
since bright quasars are black holes accreting close to the Eddington limit, and are hosted by 

' dark matter haloes with a narrow mass range of a few IO'^/i^'Mq. Therefore the clustering 

\ of bright quasars cannot be used to disentangle lightcurve models, but such a discrimination 

would become possible if the observational samples can be pushed to significantly fainter 
limits. Overall, our clustering results for the simulated quasar population agree rather well 
with observations, lending support to the conjecture that galaxy mergers could be the main 
^ ' physical process responsible for triggering black hole accretion and quasar activity. 
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1 INTRODUCTION 

At the beginning of this century, very bright quasars powered by 
black holes (BHs) with masse s of the order of lO^M© were discov- 
ered at redshifts up to z '-^ 6 j pan et al" I l200d.l200ll) . At the same 
time. X-ray observations showed that the space density of Active 
Galactic Nuclei (AGN) peaks at z ~ 2 - 3, and AGN with high X- 
ray luminosities are more common at higher red shift with respect to 
their low-luminosity counterpart |steffen et al .ll2003l: ICowie et al 



2003; Cattaneo c& Bemardi 2003; Ueda e t al.ll2003l ; lHasinger et al' 



2005) . Heckman et al. (2004), using optical data from SDSS, found 
that at low redshift only BHs with a mass < 1O^M0 are actively 
growing. Combined, these observations suggest that supermassive 
black holes (SMBHs) grow "anti-hierarchically": the more massive 
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BHs were already in place at high redshift, and since then the ac- 
cretion activity has shifted to smaller scales. 

Understanding how this evolution of BH growth relates to cos- 
mic structure formation, how BH accretion depends on the environ- 
ment, and how BHs interact with their host galaxies, have become 
central questions of cosmology that need to be answered for a full 
understanding of galaxy formation. In fact, it has become clear in 
recent years that BHs and galaxies are linked and mutually influ- 
ence each other. This co-evolution has been explored in recent years 
through analytic, semi- analytic and fully numerical approaches in 
numerous studies (e.g.. Silk &. Rees '1998'; Kauffmann & ] 
[2000 ; Merloni etal. 2004; Di Matteo et al. 2005; Cattaneo e 
2005; Croton et al. 2006; Monaco et alJl2007l ; lMalbon et : 
IMaruUi et al..2006.,.2007.) 

Starting from lSoltaij h982l) . many papers have combined the 
present-day BH mass function with the AGN luminosity density of 
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quasars at various redshifts to conclude that most of the mass in 
todays BHs must have b een accumulated during phases of bright 
AGN acti vity (see also Yu & Tremaine' 2002': El vis et aP I2OO2I : 
iMarconi et al . 2004; Merloni & Heinz 2008). The duration of these 
highly-efficient accret ion phases co uld range from a fe w 10^ yr 
jYu & Tremairi3l2002h up to lO^yr jMarconi et ai]|2004ll . values 
that strongly depend on the BH mass range considered and on the 
assumed radiation efficiency e. In fact, the pr ecise value of this 
quasar lifetime is still an open question jMartini, 2004) . Estimates 
of the duration of i ndividual accretion events using, for exa mple, 
the proximity effect dCarswell et aIlll 982':'Baitlik et al."l988^. have 
suggested lifetimes of the order of 1 Myr ( Kirkman & Tvtler 200§). 

iHaiman & HuH j200lh and iMartini & Weinberg) | |2001I) sug- 
gested to use quasar clustering to ob tain estimates of the q uasar life- 
time (see also the seminal work of ICole & Kaisei|[l989l) . The rea- 
soning behind this conjecture is simple: if quasars are strongly clus- 
tered, their hosts must be rare objects, and therefore they must also 
be long events in order to account for the total quasar luminosity 
density observed. If, on the other hand, their clustering is compa- 
rable to the clustering of small dark matter haloes, their hosts must 
be much more common, and their luminous phases must therefore 
have short duration. 

The advent of wide-fi eld surveys like SPSS and 2dF quasi- 
stellar object (2dFQSO) jYork et alj l200d : ICroom et al.l |2004|) 
with their observation of thousands of quasars has allowed a de- 
tailed investigation of th e clustering pr operties of accreting BHs. 
ICroom et al. ( 2005) and iPorciani et al.l ^2004.) calculated the cor- 
relation function of quasars observed in 2dF in the redshift range 
0-5 ^ z 2. Both groups found that the clustering strength is an 
increasing function of redshift, but that it does not depend signifi- 
cantly on quasar luminosity. The inferred values of the bias would 
suggest that quasars of the observed luminosities are hosted by 
haloes of a few lO'^/i^'M©, which remains approximately con- 
stant with redshift, since haloes of a fixed mass are progressively 
more clustered towards higher redshi ft (see also G razian et al. 
2004h. Following l.he app roach of iHaiman & Hull 1200 ll) and 
Martini & Weinber i ( l200lh . the estimated quasar lifetime would 
be a few 10^ yr, reaching ~ lO^yr at the highest observed redshifts. 



have confirmed these results 4Shen et al. 2007; Myers et al. 20071. 


Coil et al. 


20071: Ida Angela et al.ll2008l: IPadmanabhan et al.ll200a: 


Ross et al 


|2009|). However, the magnitude range covered bv these 



surveys is typically quite narrow, and this may explain the lack of 
evide nce for a significant dependence of clustering on luminosity. 
When lshen et al . (2008) analyze the clustering of the 10% brightest 
objects of their sample, they find that these quasars have a higher 
bias compared to the full sample. 

Using hydrodynamical simulations of i s olated galaxy mergers 
('Springel et al.' l2005l : iDi Matteo et alJl2005l) , lHopkins et al.l ^20051) 
studied the luminosity distributio n of accreting BHs, wh ose activity 
is triggered by the merger event. iHopkins et al.l j2005t) found that 
the luminosity distribution of the simulated AGN is equivalent to a 
highly efficient accretion phase (with very high Eddington ratios), 
followed by a decaying phase where AGN spend most of their life. 
During this extended period, they would appear as faint AGN, even 
though they may, in fact, co ntain quite massiv e BHs. 

Based on these results, iLidz et al.l i [200^ explored the depen- 
dence of quasar clustering on luminosity, using an analytic ap- 
proach to connect quasars, black hole masses and halo masses. In 
a quasar model in which the bright end of the luminosity func- 
tion is populated by BHs accreting close to their peak luminosity, 
and the faint end is mainly populated by BHs accreting at low Ed- 



dington ratio, there should be no strong dependence of clustering 
on quasar luminosity, i.e., bright and faint AGN should actually 
be the same type of objects, but seen in different stages of their 
evolution. They should therefore be hosted by dark matter haloes 
of similar masses and hence exhibit similar clustering properties. 
Assuming a relation between th e quasar B - band peak luminosity 
and the mass of the host haloes, iLidz et al.l f2006l) tested this pre- 
diction, and indeed found that only a narrow range of halo masses 
should host active quasars, with a median characteristic mass of 
Mjiaio ~ 1-3 X lO'^Mpj. As pointed out by the authors, only future 
surveys that will be able to observe the faint quasars in their qui- 
escent stage will be able to test this picture, and to rule out the 
alternative hypothesis of luminosity-dependent quasar clustering. 

In the present work we explore the properties of quasar clus- 
tering using a semi-analytic model for galaxy formation and BH 
accretion developed on the outputs of the Millennium Simulation 
ISpringel et al. 2005). Compared to other theoretical work, we do 
not have to make assumptions about the halo population hosting 
BHs nor about the relation between the halo mass and the quasar 
luminosity (or BH mass), since they are a natural outcome of the 
simulation of the galaxy formation process. However, we have to 
make assumptions about the physics of BH accretion, and what 
triggers AGN activity. In this work we are especially interested in 
testing the assumption that galaxy mergers are the primary phys- 
ical mechanism responsible for triggering accretion onto massive 
BHs. To this end we explore the simulation predictions for quasar 
clustering and the quasar luminosity function obtained with a pure 
Eddington-limited lifetime model and a mod el that includes a low- 
luminosi ty accretion mode as described bv iHopkins et 
(see also lMarulli et al . 2008, hereafter Paper I). 

After an introduction to our methodology and a review of 
some basic properties of our simulated AGN population (Section 
2), we show their clustering properties in Section 3, where we also 
compare our results with the most recent observational optical data 
available. In Section 4, we show the relation between luminous 
BHs, quiet BHs and their host haloes. Finally, we summarize and 
discuss our results in Section 5. 



2 MODELS FOR BLACK HOLE ACCRETION AND 
EMISSION 

In this Section, after a short overview of our semi-analytic model, 
we describe the different phases of BH growth and emission 
adopted in our model. We then review some basic properties of 
the simulated AGN population; further details are given in Paper I. 



2.1 From dark matter particles to galaxies 

The semi-analytic model used in thi s work is run on th e out- 
puts of the Millennium Simulation dSpringel et al.l 1 20051) . This 
is an N-body simulation which follows the cosmological evolu- 
tion of 2160"^ ~ 10'^ dark matter particles, each with mass ~ 
8.6 X lO^/i^'Mp), in a periodic box of 500/i^'Mpc on a side. 
The cosmological parameters used in the simulation are the 
ones of the WMAPI & 2dFGRS 'concordance' ACDM frame- 
work: i2,„ = 0.25, f^A = 0.75, Q% — 0.9, Hubble parameter h = 
j/p/IOO kms~'Mpc ~' = 0.73 and primordial spectral index « = 1 
dSpergel et alj20()3l) . 

The merging history trees extracted from this simulation de- 
scribe the detailed formation history of DM haloes and their sub- 
haloes, identified with a friends-of-friends (EOF) group-finder and 
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an ex tended version of the SUBFIND algorithm jSpringel et alj 
l200lh . respectively. Using the trees as basic input, our semi- 
analytic code describes the baryonic processes of galaxy formation 
and allows the prediction of galaxy properties in a large cosmolog- 
ical volume. 

The pre sent work is based on t he galaxy formation mode l 
described by ICroton eTall ( l2006h and IPe Lucia & Blaizo3 dlOOTh . 
which we extended to follow the details of BH accretion and the 
lightcurves of AGN. We refer the reader to these papers for a full 
description of the baryonic physics which describes the evolution 
of galaxies, their stars and their gas. Below, we describe only the 
prescriptions that directly relate to our study of the evolution of 
SMBHs. 



2.2 Creation and accretion of BHs 

In this subsection we discuss our modeling of the physical pro- 
cesses responsible for BH accretion. In the semi-analytic model, a 
fraction of the mass of a halo is assigned to baryons in the form 
of hot gas, which as time evolves, will cool and form a galaxy. We 
also add a 'seed' BH of very small mass to each newly formed halo. 
As galaxies evolve their central BHs are allowed to grow through 
mergers with other BHs and through gas accretion during the 'radio 
mode' and during the 'quasar mode'. The quasar mode is the phase 
during which BHs accrete most of their mass, and during which 
BHs can shine as bright AGN. We will therefore mainly concen- 
trate most of the discussion on this phase, and we will describe 
which physical process might be responsible for triggering this ac- 
tivity. 



2.2.1 BH seeding 

The origin of primordial massive BHs is still subject of intense de- 
bate: SMBHs seeds could grow out of the remnants of Pop III stars 
(e.g., lMadau & Reesll200lMHeger & Wooslevll2002h or could have 
their origin dir ectly in the collapse of a low-angular momen tum 
gas cloud (e.g.. iLoeb & Rasid[l994l : iKoushiappas et al.ll2004) . In 
the first case the progenitors of SMBHs would have a mass of the 
order of 10^ — 10-^ M©, much less than what could be the outcome 
of low-angular momentum gas collapse (lO^Mm). 

Unfortunately, due to exponential growth during accretion, it 
is very difficult to use the local population of massive holes to re- 
cover information about their original mass before the onset of ac- 
cretion. On the theoretical side, simulations are being carried out to 
investiga te which model for massive BH formation is most plausi- 
ble (e.g.. lBromm & Loebll2003l ; I Alvarez et"ai] |2008). Observation- 
ally, these models for primordial BHs will hopefully be tested in 
the near future either directly through gravitationa l wave detection 
jSesa na et al. "2005';' Koushiappas & Zentner '2006"), or indirectly by 
lookin g at the e ffect that primordial BHs might have on reionization 
(e.g.. iRicotti et alj2005l : lRipamonti et al.ll2008f) . 

As in Paper I, we assume here that every newly-formed galaxy 
hosts a central BH of 10^ Mq. This seed BH may then start ac- 
creting through the processes described below. Note however that 
a much larger seed would only influence the BH evolution in our 
model at very high redshifts, but it would not influence the results 
in the redshift range of main interest in this paper, simply because 
the large growth factor soon cancels any information about the seed 
mass. 



2.2.2 Radio mode 

When a static hot halo has formed around a galaxy, we assume that 
a fraction of the hot gas continues to accrete onto the central BH, 
causing low-level 'radio' activity in the galaxy center For clarity, 
this phase, which is called in jargon radio mode because it is as- 
sociated with the activity of radio galaxies at the centre of galaxy 
clusters jBest et al]|2005[) . does not i nclude the powerful emission 
of FRII radio loud QSOs. Following Icroton et al.l ( l2006h . the BH 
mass accretion rate during these phases of radio mode activity is 
postulated to scale as follows: 



A^BH.R = KaGN 
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0.1 ) V200kms-1 



(1) 



where Mbh is the BH mass, /j,ot is the fraction of the total halo mass 
in the form of hot gas, Vyir is the virial velocity of the halo and k^gn 
is a free parameter set equal to 7.5 x lO^^MQyr^' in order to repro- 
duce the turnover at the bright end of the galaxy luminosity func- 
tion. Since /hot is approximately constant for Vvii- It ISOkms^', 
the dependence of Mbh.r on this quantity has little effect. Note 
that the accretion rate given by equation l[T) is typically orders-of- 
magnitude below the Eddington limit. In fact, the total mass growth 
of BHs in the radio re lative to the quasar mode (discussed below) 
is negligible l lCroton et al. 2006). 

It is also assumed that radio mode feedback injects energy ef- 
ficiently into the surrounding medium, which can reduce or even 
stop the cooling flows in halo centers. The mechanical heating 
generated by this kind of BH mass accretion is parameterized as 
i'BH = e^BHf^! where e = 0.1 is the accretion efficiency and c is 
the speed of light. The heating modifies the infall rate due to cool- 
ing according to: 
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For consistency we never allow m'^^^^-^ to fall below zero. In this 
scenario, the effectiveness of radio AGN in suppressing cooling 
flows is greatest at late times and for large values of the BH mass, 
which is required to successfully reproduce the luminosities, colors 
and clustering of low-redshift bright galaxies. 

2.2.3 Quasar mode 

This is the phase during which BHs accrete cold gas and build up 
most of their final mass. This phase has recently acquired the jar- 
gon name quasar mode because it is only through the very efficient 
accretion of cold gas that a BH can shine as a bright AGN, but we 
stress that this phase is also meant to include accretion of cold gas 
at low Eddington ratios. 

The tight re lation observed locally between BH mass and the 
host bulge (e.g ., Magorrian et all 199^ ; Ferrarese & MerritllboOOl : 



iTremaine et al JI2OO2I : iMarconi & Hun j|2003l) suggests that bulges 
and BHs might form during the same events and/or they strongly 
influence each other as they evolve. Simulations have shown that 
during mergers of gas-rich disk galaxies gas is channeled toward 
the nuclei of the merging galaxies through gravitational torques 
l lBames&Hemauisllll996l) . and this process can indeed be re- 
sponsible for the formation of bulges as well as for BH accretion 

(Springel et al. 2005; Di Matteo et al. 2005jL 

Based on these results, and following iKauffmann & Haehneltl 

( l2000h . in the present work we assume that the quasar phase is trig- 
gered by galaxy mergers. In practice, during merger events, we as- 
sume that the BHs hosted by the merging galaxies instantaneously 
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Figure 1. Median accreted gas AMbh.q relative to tlie final BH mass for 
eacli accretion event, for three different final mass bins. Tlie filled contours 
enclose the 25 and 75 percentiles. 



coalesce and form a single BH whose mass is the sum of the pro- 
genitor BHs, and that this resulting BH starts accreting a fraction of 
the available cold gas. In Paper I we found that we need to feed BHs 
more efficiently at high redshifts in order to build massive BHs by 
z = 5 without invoking super-Eddington accretion or much more 
massive seed masses. In this work we assume that the amount of 
cold gas accre ted during each merger depends linearly on redshift 
( ICrotonl|2006h : 



AMbh.q = 



rl 

■/BH "^cold 



(1- 



(3) 



(280kms-i /V 

where nicoid is the total mass of cold gas in the final galaxy, Zmerg is 
the redshift of the merger and 



./inerg — /merg ('^^sat/'^^central) ; 



(4) 



where /merg ~ 0.02 is a normalization parameter chosen to match 
the observed local Mbh — Msuige relation and "Jsat/mcentral is the 
mass ratio of the merging galaxies. 

In Figure [T] we show, as a function of redshift, the median ac- 
creted mass AMbh,Q! relative to the value of the BH mass at the 
end of a single accretion event, for three final mass bins. Small- 
mass BHs accrete efficiently at all epochs (higher curve), whereas 
BHs that, at the end of the accretion event, end-up in massive ob- 
jects (lower curve) accrete most of their mass at early times: at low- 
redshifts, the amount of 'new' gas accreted is relatively small com- 
pared to the mass already acquired. This behavior is in agreement 
with the apparent 'anti-hierarchical' growth of BHs: observations 
in the soft and hard X-rays have shown that the number density 
of bright AGN declines with decreasing redshift, w hile the den- 
sity of fainter active nuclei shows the opposite trend dCowie et al.l 
20031: ISteffenet alJl2003l : lUeda et alj|2003l : iHasinger et alj|2005h . 
Heckman et al.l ( l2004h used the emission lines of type 2 AGN ob- 
served with SDSS to investigate whether the decrease of the space 
density of bright objects is simply due to a decrease in the accre- 
tion rate or a decrease in the typical mass of actively growing BHs. 
These authors found that the typical mass of BHs that are today ac- 
tively accreting is < 10^ Mq, and that larger BHs are experiencing 
little accretion. 

In Paper I we showed that, at z = 0, this model for BH ac- 
cretion i s able to reproduce not only the observed Mbh ^MBuige 
relation jHiiring & Rixl2004) , but also other scaling relations, such 
as the ones between the BH mass and the galaxy central veloc- 
ity dispersion or color jMarconi & Huntll2003l : If errarese & Ford 
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Figure 2. Differential BH mas s density at z = (re d thick line) compared 
to the observational estimate o f lShankaretaDj2004) (solid black line, with 
errors enclosed in the grey shaded area). 



I2OO5I) . The z = differential mass density of our simulated BHs 
is shown in Figure |2] compared with the observational estimate of 
Shankar et al. (2004). The corresponding local mass density (for 
our cosmology with h = 0.73) i s Pbh = 3.35 x 10^ Mr Mpc^-', 
which is in good agreement with lGraham & DriveJ ( l2007h (we re- 
fer to these authors for a summary of the values quoted in the liter- 
ature). 

To study the redshift evolution of the BH population, it is im- 
portant to not only to consider the evolution of the BH mass, but 
also to relate this to the radiation output of the accretion. If we 
are interested in the instantaneous brightness of a quasar, we not 
only need to calculate how much mass it accretes, but also how 
long this takes. In other words, we need to model the lightcurve 
of individual phases of quasar activity. In Paper I we introduced 
and tested different models for the AGN lightcurve, and we com- 
pared our results wit h the AGN bolometric luminosity function of 
iHopkins et al.l( l2007l) . We here briefly describe the lightcurve mod- 
els adopted for the present study. 

At any given time, the bolometric luminosity emitted by an 
accreting BH is given by 

Lbol(0 = eMaccr(Of^ = T— ^Bh(Oc^ 



= ./Edd(O^Edd(0 =./Edd 

ffidd 



(5) 



where e is the radiative efficiency, Ledd is the Eddington luminosity, 
,/Edd is the fraction of Eddington luminosity emitted, and rgdd = 
Qjc/ {4nmpG) ^ 0.45 Gyr (note that we are here considering only 
the luminosity emitted during the quasar mode phase, thus ignoring 
the contribution from Mbh.r)- If> at any given time, the radiative 
efficiency and the Eddington ratio are known, the accretion rate is 
given by: 



rflnMBH(f) = 



dt 



(6) 
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Figure 3. Bolometric luminosity function assuming Eddington-Iimited accretion (Mod I, blue-daslied curve), or Eddington-limited accretion followed by a 
quiescent phas e of low lumino sity (Mod II, green-solid curve), with errors calculated using Poisson statistics. The luminosity functions are compared with the 
compilation of lHopkins et alj fcOOTh (grey points with best fit given by the grey band). 
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Figure 4. Probability distribution of fEdd. as a function of BH mass and 
redshift. The limits in the BH mass bins are shown in the first panel in units 
of Mq. At high redshift, most of the BHs accrete at the Eddington limit. 
Today, only the smallest BHs are experiencing efficient accretion. 



where fef(f) 
/Edd = !)• 



is the e-folding time (fef = fsalpeter if 



l-E ./Edd(') 

For simplicity, we assumed a constant radiative efficiency 
e = . 1 (average value for a thin accretion disk, Shakura & Sunyaev 
Il973h , and we explored different models for the time-evolution of 
/Edd- In this work we choose not to explore all the four models dis- 
cussed in Paper I. Instead, we will focus on two of them, which we 
regard as representative cases. The first one illustrates the simple 



case of an AGN that shines at the Eddington luminosity. It repre- 
sents a very simple model commonly used in the literature, that we 
regard as a reference case, despite the fact that, as shown in Paper 
I, fails to reproduce the AGN luminosity function at low and high 
redshifts. The second model is very close to the model called 'best' 
in Paper I and illustrates the impact of adopting a non trivial AGN 
light-curve, motivated by numerical experiments. As discussed in 
Paper I this second model provides a better fit to the AGN luminos- 
ity function. In what follows, we present a more detailed descrip- 
tion of the two models: 

• Model I: /Edd(') = const = 1. This is the simplest case, in 
which we assume that, when active, BHs accrete and radiate at the 
Eddington limit. 

• Model II: Here we assume that BHs undergo an Eddington- 
limited phase that leads to a peak luminosity L^eak^ which is then 
followed by a long quiescent ph ase at progressively lo wer Edding- 
ton ratios. Following the work of lHopkins et ai]l l2005h . we assume 
that in this long quiescent phase the average time that an AGN 
spends in a logarithmic luminosity interval can be approximated 
by: 



d/ 



dlnLboi 



= \a\t<) 



(7) 



where = tQ{L' > IO^Lq) and to(L' > L) is the tota l AGN life- 
time above a given luminosity L. lHopkins et al found from 
merger simulations that tg ~ 10^ yr over the range IO'Lq < L(,oi < 
Lpeak; here, we assume always tg = lO^yr. In the range IO^^Lq < 
^peak ^ 10''*Lf^. [Hopkins et al.l j2005l) also found that a is a func- 
tion of only the AGN luminosity at the peak of its activity, Lpeak^ 
given by a = -0.95 + 0.321og(Lpeak/10'^LQ), with a = -0.2 as 
an upper limit. 

In this scenario, the peak luminosity Lpg^k reached at the end of 
the first accretion phase is LEdd(^BH,peak)> where 

(8) 



A^BH,peak=MBH(fin)+J ' AMbh.Q ■ ( 1 - e) 



Here MbhC'ih) is the BH mass at the beginning of the accretion. 
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AMbh.q is the fraction of cold gas mass accreted, and J sets the 
fraction of gas that is accreted during the Eddington-limited phase. 
After this first phase, the BH keeps accreting the remaining cold 
gas at a progressively slower rate, as described by equation Q. In 
Paper I we set J = 0.7, a value that balances the needs of effi- 
ciently building-up massive BHs and of explaining low-/Edd BHs 
in the local universe. Most of the available gas is therefore accreted 
during the Eddington-limited phase, and the lightcurve model in- 
troduced bv .Hopkins et aU ( i2005i) is used to describe only the qui- 
escent phase. 

A direct comparison of the luminosity functions obtained us- 
ing Mod I and Mod II is shown in Figure[3] Mod I and Mod II give 
a similar population of high-luminosity AGN: bright AGN are al- 
ways produced by BHs accreting close to the Eddington limit. At 
high redshifts, the faint-end of the luminosity function produced by 
the two models is very similar as well, suggesting that at very high 
redshifts BHs of all masses typically accrete at /^dd = 1. It is in 
the faint-end of the luminosity function at low redshifts where the 
two models predict a different behavior for the AGN luminosity: 
only Mod II (with J = 0.7) is able to fit the low-redshift faint-end 
of the luminosity function, implying that a model in which BHs 
experience long, quiescent accretion phases can indeed explain the 
number density of low-luminosity AGN at low redshift. This is be- 
cause in Mod II the average lifetime of AGN is much higher (a 
larger fraction of time is spent at low luminosities); it is therefore 
more probable to observe, at a given redshift, an AGN shining at 
low luminosities. For a more detailed discussion on this, we refer 
again the reader to Paper I. 

We have already mentioned that observations indicate that the 
more massive BHs have accreted most of their mass at early times, 
whereas in the loca l universe BHs with a mass < 10^ M0 are ac- 
creting efficiently ( Heckman et al. 2004). These resul ts hav e been 
confirmed more recentlv by [Netzer & Trakhtenbrolj ( l2007h . who 
found that at all redshifts /gdd is smaller for larger mass BHs. 
Similar compilations that use emission lines to estimate Edding- 
ton ratios have shown that the /^dd of quasars seems to be log- 
normally distributed, with a peak aroun d /^dd ~ 10^' — lO^*^'* 
jKoUmeier etall 12004 [siien et alj I2OO8I) . In Figure H we show, 
for Mod II, the redshift evolution of the probability distribution 
^(/EddlA^Bn) of the Eddington ratios, given the BH mass. At high 
redshifts all BHs accrete close to the Eddington limit. At lower red- 
shifts instead only the smaller BHs are accreting at high Eddington 
ratios, while the more massive ones accrete at much lower rates. 
Note that this figure includes all active BHs from our simulation, 
and therefore a direct comparison with observed data is not possi- 
ble. We postpone a more detailed analysis of this point to a future 
work, but we stress that a model with a quiescent phase could ac- 
count for the low-reds hift behavior of the more m assive BHs (see 
also the recent work o f lHopkins & Hernquislll2008h . 



3 CLUSTERING PROPERTIES 

In this section we discuss the clustering properties of our simulated 
AGN sample. We first compare the predicted two-point correlation 
with the autocorrelation of the DM particles. We then compare the 
AGN clustering properties with the clustering of the dark matter 
haloes of the Millennium simulation, and in particular examine the 
differences between Mod I and Mod II. We then explore the lumi- 
nosity dependence of the clustering of the global AGN population 
and of an optically-visible sub-sample. Finally, we directly com- 



pare the clustering of our simulated L* quasars with recent obser- 
vational results. 

3.1 Brief description of the correlation parameters used 

We use the standard definition of the two-point spatial correlation 
function as the excess probability for finding a pair of o bjects at a 
distance r, each in the volume element dVi and AV2 (e.g.. |Peacockl 
il999l) : 

dP = )i2[i+^(,.)]dVidV2, (9) 

where n is the average number density of the set of objects under 
consideration. 

The clustering length ro is defined as the scale at which the 
two-point correlation function is unity: ^(rg) = 1 (i.e., the scale 
at which the probability of a pair is twice the random). At scales 
between ~ I A^'Mpc up to few tens of Mpc the observed quasar 
correlation function can be approximated by a power-law, usually 
expressed as: 

^(')={£) ' ■ (10) 

To calculate rg, unless otherwise stated, we will fit the two-point 
correlation with such a power-law in the range 1 <r<20/i"'Mpc 
(see the next subsection for details on this). 

Finally, the bias between two classes of objects (e.g., AGN 
and dark matter) is defined as the square-root of the ratio of the 
corresponding two-point correlation functions: 

, _ /^agnW 

oagn.dm = v "e tt- (h) 

V ^dm('-) 

In principle, an accurate determination of the 'cosmic- 
variance' errors of these quantities as measured from the simulation 
could be calculated from the variance over many different realiza- 
tions of the universe. As we have only one simulation as large as the 
Millennium run at our disposal, this is not practical. A reasonable 
alternative is to estimate the errors by subdividing the whole Mil- 
lennium volume into sub-cubes, and then by calculating the vari- 
ance among the measurements for each of these sub-volumes, an 
approach we will follow here. 

In order to directly estimate the impact of the cosmic vari- 
ance in the predicted AGN clustering, it is necessary to model the 
AGN properties in mock samples designed to match the re al ones. 
We h ave followed this approach in a parallel work ( Maru iu"et al.l 
( l2009i) . submitted to MNRAS), where we have used the same semi- 
analytic model presented here to construct mock AGN catalogues 
mimicking the Chandra deep fields. 

3.2 AGN and darlt matter clustering 

We here show the results for the shape of the two-point correlation 
of the AGN sample, comparing it to the one of the Millennium 
dark matter particles. For simplicity, we present only the results 
obtained with Mod II, since the conclusions of this subsection are 
independent of the assumed model for the lightcurve. 

In the top panels of Figure |5] we plot the two-point correla- 
tion of the DM particles (dotted line) and the two-point correla- 
tion of faint (LboI < 10'°Le) and bright (LboI > lO^Lg) AGN 
(dashed lines), at three different redshifts. As can be seen at a 
glance, the main difference between ^dm('") and ^agn('') h^s in 
the normalization, they are substantially biased relative to each 



The clustering of quasars in semi-analytic models 7 



r [h ^Mpc] 
1 10 




10 



1 10 
r [h"^Mpc] 



10 



Figure 5. Upper panels: two-point correlation function of the Millennium dai'k matter particles (dotted line) compared with the congelation of the AGN 
population, divided into a faint and a luminous sub-sample, depending on their bolometric luminosity (as indicated in the first panel. Central panels: bias 
between the two AGN samples and the dark matter as a function of scale. Lower panels: two-point coirelation from the upper panels divided by a power-law 
fit. If ^(r) was a perfect power-law, the ratio should be constant with scale and equal to unity (dashed horizontal line). We refer to the text for a description of 
the errors. 



other. This bias ([^agn('")/^dm('')1'''^) is plotted in the next set 
of panels of Fig. |5] The bias is approximately scale-independent 
(at least in the range 1 < r < 20 /j^'Mpc), and its average value 
increases with redshift. The error CTiog^^|,j^(r) of the two-point cor- 
relation is here the variance (in log-space) of the two-point corre- 
lation functions calculated in eight sub-volumes. The errors on the 
bias have been calculated assuming negligible error for DM au- 
tocorrelation. By error propagation, the error on the bias is then 
0,,(r)= ^^(r)Oi„g^^^Jr)(lnlO)/2. 

Finally, in the lower panels of Fig. |5] we show how the two- 
point correlations deviate from a power-law, that is, we divide the 
calc ulated by the fit c alculated using eq. l IlQI l. As is well known 
(e.g. ISpringel et aLll2005l) . the DM correlation function deviates 
from a pure power-law at low and intermediate scales. The AGN 
correlation function shows a similar shape at intermediate scales 
(r ~ few h^^Mpc), but not at small scales, where the AGN two- 
point correlation function is a significantly 'better' power-law that 
of the DM. This is highly reminiscent of the findings for the clus- 
tering of galaxies dSpringel et al.ll200^ . 

As we will again see in the next subsection, the lack of a strong 



correlation signal at small scales is due to the fact that our BHs 
accrete gas and can shine as bright AGN only after merger events, 
which, in our model, happen mainly in the central galaxies of dark 
matter haloes (whose mean separation is fii 1 /i^'Mpc). Also note 
that each of our mergers lights up only one BH, the merged BH 
of the two progenitor galaxies, i.e. our model does not account for 
the possibility that the two BHs exhibit activity as a close quasar 
pair already prior to coalescence. Also, as a BH is still accreting 
cold gas, it can happen that its host halo merges with another halo 
which could have at its center another accreting BH. This is also 
why the correlation power at scales < 1 /i^'Mpc is non-zero, but 
negligible. In a forthcoming paper we will compare a pure merger- 
triggered AGN scenario, with a model in which the possible galaxy 
disk instability also could contribute in feeding BHs. In this last 
case we expect a larger AGN halo occupation distribution (number 
of AGN in a single halo), and a different behavior in the small-scale 
clustering regime. 
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Figure 6. Two-point correlation function for the AGN sample compared 
to the two-point coirelation of the Millennium FOF haloes, at vaiious red- 
shifts. The AGN are divided into 4 luminosity bins (depending on the bolo- 
metric luminosity), whereas the haloes are divided into two bins, depending 
on the value of their virial mass in units of /i^'M©. The AGN in this figure 
have been obtained using Mod II for the lightcurve. In Figure [8] the main 
difference in the coiTelation between the two models is highlighted. 



3.3 AGN and halo clustering 

In this subsection we compare the AGN clustering with the clus- 
tering of the Millennium haloes. In our model, BHs are allowed 
to accrete cold gas only during merger events, which are experi- 
enced mainly by the galaxies sitting at the centers of FOF haloes. 
As discussed above, only a small fraction of AGN can be hosted by 
satellite haloes. Due to this uncertainty in the quasar pair regime, 
we focus in the present work on the clustering on intermediate and 
large scales, and we refrain from drawing strong conclusions from 
the results at scales much smaller than the average halo separation. 

In Figure |6] we show at different redshifts the two-point cor- 
relation function of the AGN population, divided in four luminos- 
ity bins depending on their intrinsic bolometric luminosity. This is 
compared with the two-point correlation of the FOF haloes, divided 
into two bins according to their virial mass. The AGN shown in this 
figure have been obtained using Mod II for the lightcurve. The cor- 
responding correlation lengths are shown in the lower panel of Fig- 
ure [7] In the upper panel of the same figure the correlation lengths 
of the AGN obtained using Mod I are plotted, also divided in four 
luminosity bins. In the analysis of the results, we allow the expo- 
nent Y of the power-law ansatz for the correlation function to vary 
in each fit. The values of and y for the two models thus obtained 
are given in Tables [T] and |2] We also fitted the brightest bin with a 
quadratic function (ro(z) = a + fc (1 + z) + c (1 + to compactly 
summarize the results, and the values of the coefficients are given 
at the end of each table. 

Comparing the values of the correlation lengths obtained with 
the two models, we do not find significant differences, except for 
the faintest AGN (LboI < lO'^Ltri). An enlarged view of the behav- 
ior of the correlation strength of these faint objects obtained with 
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Figure 7. Correlation length as a function of redshift of the AGN sample 
divided in four bolometric luminosity bins, compared with the correlation 
length of the length of two mass bins of the Millennium FOF haloes. The 
AGN have been obtained using Mod I (upper panel) or Mod II (lower panel) 
as lightcurve models, respectively. Fits to the brightest bins are shown with 
the dotted curve. 



Mod I (solid blue curve) and with Mod II (dotted green curve) is 
shown in Figure [8] While at high redshifts there is hardly any dif- 
ference between the two models, at low redshift the faint objects 
obtained with Mod II are much more strongly clustered. This is 
because most of the population is composed of large BHs that are 
accreting at low /gdj (as shown in Figure|4j and that are hosted by 
large haloes. In the lower panel of Figure [8] we see that the corre- 
lation length of the faint objects obtained with Mod II is compara- 
ble to the ones of haloes with Myir ~ lO'^ — lO'^M.j, while faint 
objects obtained with a pure Eddington-limited accretion model 
are sitting in haloes of much lower mass. Observational cluster- 
ing measurements have been used in recent years to estimate the 

gpical halo masses that host qu asars (e.g.. IPorcian i et al. 200^; 
-azian et al.ll2004 ICroom et al.ll2005h . This is usually done by 
comparing the bias of obser ved quasars with the halo bias obtained 
from analytical estimates jMo & Whi^ 1 19961 : ISheth & TormeiJ 
1 19991 e.g.,). In the present work, the host halo mass is an output 
of the simulation, and therefore we can directly examine the rela- 
tion between black hole mass, quasar luminosity and halo mass. In 
section|4] we exploit this for a direct study of the dark environment 
of luminous BHs. 

Based on Figure [T] it seems that the redshift-evolution of the 
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Table 1. Values of the correlation lengths shown in the upper panel of Figure [T] We also added the values of the corresponding power-law slope y. Li 
corresponds to the brighter bin, L4 to the faintest. We also give the values of the parameters of the quadratic fit done on ro for the brightest bin. 
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Table 2. Same as the previous table, this time for the AGN obtained with Mod II (lower panel of Figure|7j. 



clustering of quasars is consistent with the redshift-evolution of the 
clustering of dark matter haloes (quasars of a given luminosity re- 
side at all times in haloes of a fixed mass). Again, the only sub- 
stantial difference to this trend is for the faint objects obtained with 
Mod II: since their clustering is more constant with redshift, it im- 
plies that their typical host halo mass changes with redshift. 



3.4 Luminosity dependence of AGN clustering and 
comparison with observational data 

In this subsection we first examine the dependence of AGN cluster- 
ing on luminosity, looking at the global population, and then con- 
sidering a subsample that would be observable in the optical band. 

Observationally, quasar clustering seems not to depen d sig- 
nificantly on luminos ity jPorciani et al. 2004; Groo m et al.ii2005l : 
Ida Angela et"ai]|2008l , e.g.,). Onlv lshen et al.. qOOm found indica- 
tions of a luminosity-dependence of the clustering when they com- 
pared the two-point correlation of their 10% brightest objects with 
the rest of the sample. Figure |7] provides information on how the 
correlation length evolves with luminosity in our models. Except 
for the faintest bin (see Figure [§}, there is not a substantial differ- 
ence between the two models, as pointed out before. In both models 
we see some moderate evolution with luminosity, and in particular, 



in both cases the brightest quasar bin is substantially more strongly 
clustered than the lower luminosities. 

Note that in this analysis a very large range in luminosities 
is covered (~ 5dex in luminosity, corresponding to ~ 12.5 abso- 
lute magnitudes). Observationally, the accessible luminosity range 
is always much smaller than that. To give predictions that can be 
compared with future observations, we now extract from the global 
AGN population sub-samples of optically visible bright AGN. First 
of all, to account for obscuration, we calculate the fraction of ob- 
jects that would be visible in the optical using the 'observable frac- 
tion' from lHopkins et aL I( l2007h . This gives, as a function of lumi- 
nosity, the probability for an object to be seen in a given band: 



m=f46 



(12) 



lO^^ergs-i, 

where /4f, = 0.260 and (3 = 0.082) for the B-band. 

To convert from bolometric luminosity to B- band luminos- 
ity, we used the bolometric corrections again from iHopkins et all 
( I2OO7I) : 

^1 / , \k, 



^band 



= Cl 



^bol 

10'"L,: 



^bol 
IOIOL,; 



(13) 



where {cl,kl,c2,k2) are respectively (6.25,-0.37,9.0,-0.012) 
for the B-band. 
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Figure 8. We compare here the correlation function of faint AGN (LboI < 
IO'^Lq) obtained using Mod I (solid blue line) and Mod II (dotted green 
line). We show the result at very high redshift, where there is no difference 
in the two models, and at low redshift, were the difference becomes signif- 
icant. In the lower panel the corresponding correlation function is shown 
as a function of redshift, and the correlation of FOF haloes is shown for 
reference. 
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Figure 9. Space density as a function of redshift for four subsamples se- 
lected with B-band magnitude cuts as indicated on the plot. The solid lines 
are give the space density when the possible obscuration is taken into ac- 
count. If we allow all our objects to be optically visible, we obtain the space 
densities described by the dotted curves. The dashed line marks the point 
below which we have less than 500 objects remaining in in the Millen- 
nium simulation volume. The open diamonds are the observed values from 
iPorc iani et al. (2004), obtained in different magnitude ranges depending on 
the redshift (see text for details). The number densities obtained with our 
model using the same magnitude ranges and accounting for obscuration are 
indicated with the filled circles. 




Z 



Figure 10. Correlation length (top panel) and bias (lower panel) for the 
AGN selected using the cuts of Figure |9] (neglecting the effects of obscura- 
tion). Due to lack of enough objects, the clustering properties of the two 
brightest bins are calculated only down to z = 1.5. Our predict ions are 
plotted together with observational data (for the lShen et al T j2008l) . we m- 
cluded their lower es t imates ). For the bias, the dotted line is the predic- 
tion of Hopkins et and the short-dashed line is the best fit from 
.Croom et al.. C2005i) . 



In Figure |9] we show as a function of redshift the num- 
ber density of our simulated AGN for different luminosity cuts 
(solid lines). In order to directly compare our number densities 
with the values inferred from observational data used for cluster- 
ing measurements, we calculated in the same figure the number 
density of objects in the magnitude ranges given by Porcia ni et al.l 
{2004) at three different redshifts: the values of M^j^ and Mn,ax 
are [-25.32,-21.72] at z ~ 1.0, [-25.97,-22.80] at z ~ 1.5, and 
finally [-26.44,-23.37] at z ~ 2.0 (see their Table 1). Note that 
their value are in bj, and to convert from B to the iy-band we 
used the relation given by these authors in their Appendix 1, where 
Mb = Mj,j + 0.07. In the Figure, our po i nts are the black dots, while 
the numbers quoted by IPorciani et ( I2OO4I) are shown with dia- 
monds (the errors quoted by these authors are comparable to the 
size of the symbol, and therefore are omitted). The agreement is 
quite good, even though we slightly underestimate the number of 
bright quasars at z = 2, as expected (see the bright-end of the lu- 
minosity function at this redshift in Figure[3]). In Figure[9]we also 
show the number density of our simulated AGN for the same lumi- 
nosity cuts, but without accounting for obscuration (dotted lines). 
As described above, we account for obscuration by calculating 
for each object its probability of being optically visible and then 
by randomly extracting objects that satisfy the imposed condition. 
Since this probability is a weak function of luminosity, and since 
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Figure 11. Correlation length (top panel) and bias (lower panel) for L, 
quasars. The gray line is our prediction (with errors enclosed in the grey 
area). The observational data are the same of Figure [To] A fit to our pre- 
dicted bias as a function of redshift is given in equation )14t . 



clustering analysis is independent of random sampling, for our 
study we ignore the effect of obscuration. This allows us to push 
the analysis to brighter magnitude cuts, since for a statistically- 
accurate clustering analysis we need at least few hundred objects 
(the dashed horizontal line shows the point at which, in the full 
simulation volume, we cannot expect more than 500 objects). 

The correlation lengths of the AGN selected with these lumi- 
nosity cuts are shown in Figure[TO] We see that at low and interme- 
diate redshifts the correlation length and the bias depend weakly on 
luminosity when a narrow range of luminosities is examined. Since 
bright quasars are always powered by BHs accreting close to the 
Eddington limit, it seems difficult to use quasar clustering obser- 
vations to disentangle between different light-curve models, unless 
much larger luminosity ranges are probed. The present observations 
indicate however that, over the range of luminosities observed, 
quasars reside in haloes of similar masses. Based on our results, 
we conclude that the lack of a significant dependence of clustering 
on luminosity is not primarily a result of invoking lightcurve mod- 
els with a wide distribution of Eddington ratios, but rather arises 
because in a merger-driven scenario there is a small scatter in the 
typical halo mass hosting quasars close to their peak luminosity. 

In Figure[TO]we added observational data from several works, 
to qualitatively compare our results with observations. We stress 
though that the error bars in these figures are calculated to de- 
scribe the effect of cosmic variance as described in Section [3T1 
since we are here ignoring the effect of obscuration, thus improv- 



ing our statistic, a direct comparison with the error bars given by 
observational works is not possible. 

Most of the observe d quasars have a typical magnitude around 
jCroom et alj|200^ . with faint limits that strongly depend on 
redshift (at very low redshifts surveys can reach fainter magni- 
tudes, whereas at very high redshifts the limiting magnitudes can 
be higher than M*). At z < 1 the faintest observed magnitudes are 
Mb ~ —22, going up to ~ —24 at z ~ 2 — 3. Since each observa- 
tional study uses different magnitude cuts, we can not do a detailed 
comparison with all the observations available, but overall our re- 
sults for the values of the correlation length and the bias and their 
evolution with redshift are in good agreement with the observa- 
tional results. 

We also compared observational data with simulated quasar s 
around L», calculated using equation 9 from lHopkins et al ]( l2007h . 
and selecting objects with an intrinsic luminosity larger than 
L*/0.5dex (which corresponds to a minimum luminosity approx- 
imately 1.2 mag fainter than M^). Our predictions for the correla- 
tion length and the bias for objects as a function of redshifts are 
shown in Figure[TT] again to gether with the ava ilable observational 
data. The discrepancies with I Shen et alj ( |2008|) for the correlation 
length can be due in differences in the calculation of this quantity 
(as already mentioned, here we do n ot fix the value of y) - For the 
bias, we sho w also the best fit fro m ICroom et alj ( |2005|) and the 
prediction of iHopkins et al. I l l2007l) . The latter was probably fitted 
only up to z = 3, thus explaining the turn-over at redshifts above 
3 that seems to not be consistent with the trend shown by the ob- 
servations. A good approximation to our prediction for the bias is 
given by the fitting function 



fe(z) =0.42 + 0.04(1 +z)+ 0.25(1 +z) 



(14) 



Quasars with luminosities around L* are typically objects very 
close to their peak luminosity, therefore correspond to objects ac- 
creting at high Eddington ratios. As mentioned before, we can- 
not use these results as a sensitive test of our lightcurve models. 
However, the good agreement with observations indicates that our 
merger-triggered BH accretion model predicts a spatial distribution 
of quasars that is consistent with observations. This is a prediction 
of a consistent model of the joint evolution of dark matter, galax- 
ies and black holes, evolving ACDM initial conditions from high 
redshift to the present. While the parameters of the semi-analytic 
model had been tuned to fit the bulk z = properties of the BH pop- 
ulation and the AGN luminosity function as a function of redshift, 
information on clustering had not been considered in the construc- 
tion of the model, and therefore must be regarded as genuine model 
predictions. 



4 BHS, QUASARS AND THEIR DARK ENVIRONMENT 

In this section we explore directly the connection between BHs, 
quasars and their dark matter environment. As in our simulations 
the dark matter halo merger trees are the backbone upon which 
the baryonic component is treated, we can also use them to study 
the dark environment of our AGN. This in particular allows tests 
of the validity of the approach typically adopted in the interpreta- 
tion of observational quasar clustering results (e.g. IPorciani et al.l 
^2004; Crooniet aL 2 005,) , where the observed quasar bias is com- 
pared with the halo bias predicted by ana lytical halo models (e.g. 
IMo & Whitell996l : ISheth & Tormenll 19991) . 

The mass distribution of the haloes hosting AGN of given lu- 
minosities, /'(MHaiol'-AGN)> is shown in Figure[T2] The AGN are 
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Figure 12. Distribution of dark matter halo masses hosting faint-AGN, blight AGN and L, quasars. The vertical dashed Hne indicates the median of the 
distribution for each luminosity bin, and we refer the reader to the legend on the plot for details in the color/pattem-coding. The AGN have been obtained 
using Mod I (upper panel) of Mod II (lower panel) for the lightcurve, respectively. 



here sub-divided into a faint and a bright sub-sample, depending 
on their bolometric intrinsic luminosity. The cut in bolometric lu- 
minosity is here L*, calculated in the same way as for section [J!4l 
Based on the results on the Eddington ratio distribution (see Figure 
|4]l and on the clustering, we expect the distribution of the masses 
of the haloes hosting bright AGN to be similar both for Mod I and 
Mod II. The main difference should be in the distribution of haloes 
hosting faint AGN: in the Eddington-limited model, the faint AGN 
population is composed of small-mass BHs accreting at Eddington, 
whereas in the model that includes a long quiescent phase the faint- 
AGN population at low redshifts includes also quite massive BHs 
accreting at low Eddington ratios. 

In Figure [12] we indeed see that for Mod I there is a direct 
proportionality between the luminosity of the AGN and the mass 
of the host halo: the brighter the AGN, the larger the BH and the 
host halo. Instead, for Mod II most of the low-luminosity AGN at 
low redshifts are hosted by more massive haloes, i.e., massive BHs 
accreting at low Eddington ratio. In the same figures we also plot 
the mass distribution of haloes hosting L, quasars. To get a large 
enough sample, at any given redshift we included objects in a range 
of ±0.5 dex around L,. The similar behavior of haloes hosting 
quasars in both models suggests that L* objects are mainly BHs 
accreting close to the Eddington limit. 

The mean values of the distributions are shown as a function 
of redshift in Figure[T3] In recent years many groups have analyzed 



the clustering properties of quasars to estimate the typical mass of 
their host haloes, at l ow- jPadmanabha n et al. "2008"), intermediate- 



Croometal. 2d05l: IPorciani et al. 2004; da Angela et al. 200i; 



( 

iMverset all 1200% " and high- jShenTt al.l I2007i) redshifts. These 
works used quasars observed with SDSS and 2dF, with a typical 
luminosity around L» (except for the very high-redshifts measure- 
ments). The masses of the dark matter haloes hosting quasars es- 
timated by these groups are overplotted in Figure [13] Almost all 
these estimates are in the range predicted by our model: the typical 
halo mass hosting L, quasars seems to grow up to z ~ 1.5 — 2, and 
then it decreases again at higher redshifts. To compactly represent 
our simulation results, we fitted our results with a cubic function 



^halo = "0 +«1Z + «2Z +a^t 



(15) 



with ai = [11.873;0.944;-0.318;0.026] for the second panel of 
Figure [13] (the values of these coefficients are similar for the fit 
of the curve of the upper panel, which we omit for brevity). 

Our results for bright qu asars (objects arou nd L^ f) are also con- 
sistent with the estimates of iLidz et al.l ( l2006b and iHopkins et al.l 
( |2007() , who calculate that the typical mass of haloes hosting 
quasars is ~ 4 x 1O'^/i^'M0. These authors argue that bright and 
faint quasars are the same type of objects but seen in different evo- 
lutionary states, and therefore their typical host halo mass should 
be similar. Since only the brightest quasars are objects accreting at 
high /Eddi only for these objects we expect a tight relation between 
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Figure 14. Relation of LboI of the AGN versus dark matter halo mass. In 
the upper panel, BHs accrete according to the Mod I lightcurve, while in the 
lower panel the predictions are produced using Mod II. While all very bright 
objects are BHs accreting close to the Eddington limit, the main difference 
between the two models Ues in the faint objects, where we have a dense 
population of faint AGN hosted by large haloes (the light-green open circles 
in the lower panel refer to AGN in the quiescent phase). For reference, the 
dashed line marks the Eddington luminosity con'esponding to a BH mass of 
lO^'Mm 



Figure 13. In these two panels we show the redshift evolution of the median 
mass of dark matter haloes hosting AGN of different luminosities from the 
previous Figure ii2\ . For clarity in the plot, we only show the values ob- 
tained for objects with Lgoi < L, and with Lgoi ~ L,. The dotted black curve 
shows the best fit to the evolution of the typical host mass of L, quasars. 
The contours indicate the 25 and 75 percentiles. We overplot here estimates 
obtained by different groups who examined the clustering properties of ob- 
served quasars (see legend on the plots). 



the instantaneous luminosity and the host halo mass. The relation 
between AGN luminosity and halo mass is shown in Figure [T4l 
Indeed, only for the very bright quasars there is a direct proportion- 
ality between luminosity and halo mass. These are objects that are 
close to their peak luminosity, have accreted most of the gas avail- 
able, and at this point their BH is tightly correlated with the mass 
of the host halo (see also next figure). During the rising phase of 
the lightcurve (even if BHs are accreting at Eddington), BHs are 
not yet strongly correlated with the host halo, reflected in a lack 
of correlation between quasar luminosity and halo mass. During 
the decaying phase. Mod II produces a dense population of faint 
obje cts sittin g in massive haloes (see open circles in Figure [74l(. 

IWhite et^afj | |2008|) claimed that the very high bias observed 
for high-redshift quasars implies a small dispersion in the above re- 
lation. Es timates of high Edding t on ratios for brig ht objects at high 
redshifts l lKollmeier et al]|2006| - IShen et alJllOoj) indeed seem to 
support that for very bright object s a tight relation exists between 
quasar luminosity and halo mass jFineetal.|[2003) . However, we 
would like to point out that just looking at the bright quasar popu- 
lation it is not sufficient to distinguish between different lightcurve 
models. 



The observed scaling relations between BH masses and dif- 
ferent properties of the host galaxy have suggested the possibility 
of a more fundamental connection between the mass of the BH and 
the host system. Using measurements of stellar velocity dispersions 
and assuming a relation between this quantity and the circular ve- 
locity of th e galaxy and the BH mass. lFerrares 3( l2002l) .l Baes et al.l 
( 2003]) and lShankar etal] ( l2006h estimated how the BH mass could 
be connected to the dark halo mass in the local universe. At higher 
redshifts these estimates are of course more problematic, because 
studies of the stellar kinematics are unavailable and we also are 
not certain yet how the M^h — o relation evolves with redshift. 
iFine et al. I [2006 ) explored the relation between BHs and quasar 
host haloes at z = 0.5 — 2.5 using BH virial masses estimates from 
the width of broad emi ssion lines and DM h alo mass obtained from 
quasar clustering from lCroom et alj ( |2005|) . In Figure [15] we plot 
the Mbh — Mnaio relation for our simulated BHs. We include here 
only BHs residing in central galaxies of EOF haloes. This is be- 
cause in our model only central galaxies can merge, and therefore 
it is mainly BHs hosted by EOF haloes that can grow (the results 
of lLi et al.l ( [200^ indicate that this could be supported by observa- 
tions) . Indeed, we find a well-defined relation which gets tighter 
with decreasing redshift. In Paper I we already showed this relation 
at redshift z = and we found good agreement with o ther works 
l lFerraresell2003 : iBaes et al.l 120031: Ishankar et alj|2006l). Here we 
overplot the results of iFerrares'^ 20021) and IShankar et al] booel) 
at z = 0. 1 for referen ce; at z = 1 we overplot the zero-point in the 
relation estimated bv lFine et al] fco06) (Mbh = 1O^ '*='=°-^M0 for a 
halo of Mhaio = IO'^'^Mq) and at z = 1 and z = 2 the resul t s from 
direct hydrodynamical simulations of IColberg & di Matted ( 1200 Sh 
(for z = 2 we used their result at z = 3). 

Note that the fact that BHs need to accrete most of the 




Figure 15. Mbh — ^Halo relation for BHs sitting in central galaxies. The points are our simulated objects, an d the red line is t he best-fit assuming a linear 
relation. The filled region encloses the 25 and 75 percentiles. For reference, we show at z = . 1 the result that lFerraresd i2002l) obtained at z = assuming 
Vvii = (dash ed line), Vc = l .Si'yir (dot-dashed line) and the prescription from lBullock et all jlOOlh for this relatio n (solid line). At z = 0.1 we show also 
the result from lShankar et al] j2006l) (dotted curve). The point at z = 1 is the zero-point of this relation obtained by iFine et al J i200d) . The dashed lines at 
2=1 and at z = 2 are from lColberg & di Matted i2008l) (for z = 2 we used their result at z = 3). The horizontal dashed line marks Mbh = 10* M ;,, which 
is approximately our resolution. This plot was obtained assuming Mod I for the lightcurve, but the result does not change using Mod II, since the final BH 
masses are the same. The diamonds at z = 5 show the relation between BH mass and halo mass if BHs accreted the available mass instantaneously. 



available gas before they 'sit' on the above relation could be in- 
fluenced at high redshifts by the resolution limit of the Millen- 
nium simulation, which does not resolve low-mass haloes below 
~ 10^^ h~^MQ. We will explore this high-redshift behavior in more 
details in future work. 



4.1 Duty cycle 

The time BHs spend shining as quasars is still an open question 
(see review bv lMartini|[2004 . The definition itself of a 'quasar 
lifetime' is somewhat ambiguous. Observationally it is defined as 
the time BHs spend shining at luminosities higher than some limit 
(for quasars, the usual definition is the time an active nuclei shines 
with a B-band magnitude Mb < —23 mag). Theoretically, it can be 
defined in a simpler way as the total time a BH shines at high 
Eddington ratio. The quasar lifetime is often also simply defined 
through the duty cycle, which is given by ratio of the quasar num- 
ber density and the number density of the haloes that can host them: 
tHuhhung/nHaif , (e.g.. , Adelberger & Steidel 2005.) 



iHaiman & Hiif ( 1200 ih and iMartini & Weinberg! JioOlh sug- 
gested to use the observed quasar clustering to estimate the quasar 
lifetime, upon the assumption that a monotonic rel ation exists be- 
tween quasar luminosity an d halo mass (see also iHaehnelt et al.l 
ll998l) . lAdelberger & Steidell ( l2005h pointed out that the theoretical 
estimate of the duty cycle through clustering analysis depends on 
the Eddington ratio distribution, on obscuration and on the scatter 
in the realtion between quasar luminosity and halo mass. As we 
have seen, the assumption of a tight relation between luminosity 
and halo mass is overly simplistic for realistic lifetime models, and 
it is therefore interesting to use our simulations directly to examine 
the distribution of quasar lifetimes. 

In Figure[T6]we show the fraction of active haloes (or duty cy- 
cle), as a function of quasar luminosity, redshift and halo mass, for 
both Mod I (left panels) and Mod II (right panels). At high redshifts 
massive haloes have a very high duty cycle, i.e., most of haloes host 
a bright quasar. As expected, the duty cycle evolves more strongly 
with redshift for the more luminous AGN: by redshift z = 0. 1 only 
~ 0.1% of the more massive haloes host a quasar, and this result 
is independent on the lightcurve model assumed. Again, the differ- 
ence in the two models is in the faint AGN population: the duty 
cycle of faint objects evolves strongly with redshift and mass for 



Mod I, since at low redshift only the smallest haloes host an active 
BH. On the other side, if the AGN lightcurve includes a long low- 
level phase, then at low redshift also massive haloes are hosting a 
low-luminosity object. 

Estimates of the quasar lifetime obtained from quasar clus- 
tering suggest timescales of the order of 10^ — lO^yr , depe nding 
on the redshift. At high redshifts (z ^ 3.5). IShen et alj ( |2007[) esti- 
mated lifetimes of the order of 30 ~ 600 Myr, while at 2.9 ^ z < 3.5 
the estimated range decreases to 4 ~ 50Mvr. |Porciani et al.l ( l2004h 
suggest tq 10^ yr at z ~ 1, and values approaching lO^yr at higher 
redshifts. As we approach low red shifts and the local univers e, the 
quasar lifetimes seem to decrease: IPadmanabhan et al] ( 1200 Sl) sug- 
gest values < 10^ yr for their sample of quasars at 0.2 < z < 0.6. 
As we have shown in Figure [16] a strong evolution of the quasar 
lifetime is also expected from our models: at intermediate-high red- 
shifts our results are compatible with lifetimes of a few lO^yr, but 
the detailed evolution of the duty cycle also depends stongly on the 
range of host halo mass considered. 



5 CONCLUSIONS 

In this series of papers we investigate semi-analytic models for BH 
accretion and quasar emission in the context of a comprehensive 
galaxy formation model developed for the Millennium Simulation. 
The physical scenario for BH growth we study is ba s ed on the 
model for BH accretion from iKauffmann & Haehneltl ( l200(]h . as 
revised bv lCrotonetai] ( |2006|) . which assumes that galaxy merg- 
ers are the primary physical mec hanism responsible for efficiently 
feeding central BHs. In Paper I dMaruUi et al.ll200^ we used the 
most recent observations of the local BH population to test basic 
predictions of the model for the local BH demographics, testing 
also different theoretical models for the quasar lifetime with goal 
to reproduce the observed quasar luminosity function. We found an 
overall good agreement between the predicted and the observed BH 
properties, and that the faint-end of the observed luminosity func- 
tion can be better reproduced when a quasar lightcurve model is 
adopted that includes long quiescent accretion after an Eddington- 
limited accretion phase. 

In the present work we used the spatial distribution of active 
BHs as a further test of our model for BH accretion. Throughout 



The clustering of quasars in semi-analytic models 15 




10.0 10.5 11.0 11.5 12.0 12.5 13.0 10.0 10.5 11.0 11.5 12.0 12.5 13.0 

log M„„ [h- Mj log M„.„ [h"' Mj 



Figure 16. Fraction of active lialoes (or duty cycle), as a function of redshift, halo mass and AGN luminosity. We compare the results obtained for Mod I (left 
panels) and Mod II (right panels). 



the paper, we compared the results obtained adopting two different 
theoretical models for the quasar lifetime: pure Eddington-hmited 
accretion (Mod I), and a model in which Eddington-limited accre- 
tion is fol lowed by a l ong, weak accretion phase (Mod II), mod- 
eled after Hopkins et alj fcoOS) . The main difference between the 
predictions of the two models is in the faint-end of the luminos- 
ity function. The long low-luminosity accretion phase allowed by 
Mod II gives rise to a large population of massive BHs that at low 
redshifts are accreting at low Eddington ra tios, in agreement with 
the observational results of, for example, iHeckman et al.l ( |2004|) 
and lNetzer & Trakhtenbroll ( l2007h . who found that in the local uni- 
verse only BHs with mass < IO^Mq are expe riencing high-efficient 
accretion. As also recently pointed out by 'Hopkins & Herng uisll 
1 I2OO8) . it is only by studying the properties of the faint AGN pop- 
ulation that the quiescent phase described bv lHopkins et al. I( l2005l) 
can be tested. 

Independent of the model adopted for the lightcurve, the two- 
point correlation function of our simulated AGN can be approxi- 
mated by a single power-law in the range 0.5 '^r <2Q /i^'Mpc. 
The bias between AGN and the dark matter is a strong function of 
redshift, but, at a given epoch, it is approximately constant in the 
range 1.0 < r < 20 /i^'Mpc. As expected, the correlation lengths 
of AGN obtained with Mod I or Mod II differ only for the faint 
population: the correlation length of faint AGN obtained with Mod 
II is consistent with the correlation length of lO'^ — 1O'"'/i^'M0 
haloes, whereas faint AGN obtained with Mod I exhibit the same 
clustering as lO" - IO'^/j^'Mq haloes. 

Recent results from optical quasar surveys like SDSS and 
2dFQSO have not found ev idence for a strong dependence 
of clustering on lumin o sity jPorciani et al] | 2004l : ICroom et al.l 
20051: iMvers et al.ll2007l ; [da Angela et alj|2008l e.g..). except for 
ShenetalJ ( l200 8l) who detect an excess of clustering for their 10% 
brightest quasars. Our results are consistent with these observations 
if we consider only quasars with an intrinsic luminosity within the 
range probed by these surveys. However, if we compare the clus- 
tering properties of AGN over a very extended range of luminosity, 
then the correlation length becomes a moderately strong function of 



luminosity and the value of the correlation length of the faint pop- 
ulation in particularly is seen to depend on the lightcurve model 
assumed. The fact that the clustering of the observed quasars does 
not depend on luminosity could be explained in two ways: quasars 
of different luminosities are powered by BHs of the same mass that 
are in different stages of their evolution, and/or the typical mass 
of haloes hosting quasars is approximately constant for the lumi- 
nosity range probed by observations. From our results the second 
hypothesis seems to be clearly favoured. The mass range of haloes 
hosting quasars is narrow enough that a significant luminosity 
dependence of clustering cannot be detected with the current ob- 
servational samples, independent of the lightcurve model. 

We also directly compared the clustering of our quasars 
with the most recent observational data, and found very good agree- 
ment. Since quasars at these luminosities are objects very close to 
their peak luminosity, and therefore correspond to objects accreting 
at high Eddington ratios, we cannot, however, use these results as 
a sensitive test of our lightcurve models. Nevertheless, the good 
agreement with observations indicates that our merger-triggered 
BH accretion model predicts a spatial distribution of quasar that 
is consistent with observations. This non-trivial outcome can be 
viewed as a further success of the hierarchical galaxy formation 
paradigm, and the cold dark matter hypothesis. 

We note that a similar result for the luminosity depen dence 
of AGN clustering has been found in iMarulli et al.l | |2009|) , who 
analyzed mock AGN Chandra catalogues constructed with the 
same semi-analytic model adopted in this work. Furthermore, 
iThacker et alj ( |2009|) have recently found very similar results mod- 
eling the AGN spatial properties in an hydrodynamical simulation. 

In future work we will compare the merger-triggered quasar 
model with alternative suggestions for the physical triggering 
mechanism of quasar activity, such as disk-instabilities occuring 
in isolated galaxies. We expect that quasar clustering statistics can 
here be a potentially powerful discriminant to further constrain 
the viable physical models for the evolution of supermassive black 
holes, and their co-evolution with galaxies. 
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